03: Reporting Linear Models with Quarto

Overview

This tutorial covers how to run, inspect, and report a linear model in R. For the report portion, we will cover some key features of dynamic reporting in Quarto and how to write and render professionally formatted documents.

The Linear Model

In this section, we will walk through the process of fitting, comparing, and reporting hierarchical linear models in R. This is not a statistics tutorial, so there will be minimal detail about how to understand or interpret the output of these commands. Instead, refer to Prof Andy Field’s discovr tutorials, which are the primary teaching materials for the same content in UG teaching. All of the code and interpretation in the following section is from discovr_08, the GLM.

There are two key goals for this linear model walkthrough:

  • Create a detailed “cheatsheet” for a linear model analysis for future reference
  • Get familiar with the discovr tutorials

Of the two, the first is more important. I’d strongly recommend you open the relevant discovr tutorial and skim through it as you go so you’re familiar with what it contains. However, the following sections of this tutorial will present the same code and information with very abbreviated text, to serve as a quick reference.

You can also have them both open at once and refer to each!

Prof Andy Field’s discovr tutorials provide detailed walkthroughs of both the R code and the statistical concepts of a variety of statistical analyses. They are a good place to look first to understand what your UG supervisees or advisees have been taught on a particular topic.

The tutorials are built in {learnr}, an interactive platform for learning and running R code. So, unlike the tutorial you’re currently reading, they must be run inside an R session. We have already installed all of the tutorials in your Posit Cloud workspace.

To open a tutorial, open any project and click on the Tutorial tab in the Environment pane. You can run any tutorial from here, but the relevant one for the linear modelling we are working on now is discovr_08, “the GLM”. Scroll down to this tutorial and click the “Start Tutorial ▶️” button to load the tutorial.

Because discovr tutorials run within R, you don’t need to use any external documents; you can write and run R code within the tutorial itself. However, I strongly recommend that whenever you work with these tutorials, you write and run your code in a separate document, otherwise you will have no record of the code and output.

Data and Codebook

Today’s data is a truncated version of the TeachingRatings dataset from the {AER} package. You can load the codebook in the Help viewer by running the following in the Console:

?AER::TeachingRatings

One Predictor

Now that we have some data, we can fit our first model with a single predictor. We will do this with the very hardworking lm() function in R, standing for “linear model”.

The lm() function has a lot of options (as you might expect, given the versatility and ubiquity of linear models!), but its basic format to fit a linear model with a single predictor is very simple:

lm(outcome ~ predictor, data = our_dataset)

In this function, outcome ~ predictor is a formula expressing a (simplified version of) the linear model equation. Here, outcome is the name of the variable in our_dataset that contains the outcome or dependent variable y, and predictor is the name of the variable that contains the predictor or independent variable x.

Exercise

Create a new code chunk and use the Codebook and the lm() function to fit a linear model predicting teaching evaluation score from beauty ratings. Save the resulting model in a new object called eval_lm.

eval_lm <- lm(eval ~ beauty, data = teach_tib)

That’s it!

Model Information

Now we have a new object that contains all the information about our model. We could simply call the name of this object, but the output doesn’t tell us anything besides the actual value of the b estimates. Instead, we’ll use some useful functions from the {tidyverse} package {broom} to get the information we need.

Model Fit

To get some common measures of model fit, we can use the function broom::glance(). The output includes \(R^2\), adjusted \(R^2\), and F and accompanying statistics in comparison to the null model (the mean of the outcome alone).

Exercise

In a new code chunk, put the eval_lm object into broom::glance() to get model fit statistics.

broom::glance(eval_lm)
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic   p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>     <dbl> <dbl>  <dbl> <dbl> <dbl>
1    0.0357        0.0336 0.545      17.1 0.0000425     1  -375.  757.  769.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Helpfully, broom::glance() (and many other {tidyverse} functions) output tibbles, which means we can work with them as we already know how to do. In Tutorials 4 and 5, we’ll also learn more about changing and filtering tibbles that will make this very more useful. For now, we can simply note the information we get out of this function for future use.

Model Parameters

To get information about the b estimates for individual predictors, we can use the function broom::tidy(). We could run this function without any other information, as we did with glance() above, but we’ll change one argument here in order to get 95% confidence intervals for b in the output as well.

Exercise

In a new code chunk, put the eval_lm object into broom::tidy(). Use the argument conf.int = TRUE to obtain confidence intervals.

broom::tidy(eval_lm, conf.int = TRUE)
# A tibble: 2 × 7
  term        estimate std.error statistic   p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
1 (Intercept)    4.00     0.0253    158.   0           3.95       4.05 
2 beauty         0.133    0.0322      4.13 0.0000425   0.0698     0.196

Hierarchial Models

Next, we may want to test the addition of further predictors in the model. We can then compare the more complex multi-predictor model to the single-predictor model.

Exercise

In a new code chunk, fit a linear model with teaching evaluations as the outcome, and both beauty ratings and gender as predictors. Save the model output in a new object called eval_full_lm.

Then, obtain model fit statistics and parameters as before.

Hint: to add a new predictor, you will need to literally add it (+) to the formula.

eval_full_lm <- lm(eval ~ beauty + gender, data = teach_tib)

broom::glance(eval_full_lm)
# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic     p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>       <dbl> <dbl>  <dbl> <dbl> <dbl>
1    0.0663        0.0622 0.537      16.3 0.000000141     2  -368.  744.  760.
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>
broom::tidy(eval_full_lm)
# A tibble: 3 × 5
  term         estimate std.error statistic    p.value
  <chr>           <dbl>     <dbl>     <dbl>      <dbl>
1 (Intercept)     4.08     0.0329    124.   0         
2 beauty          0.149    0.0320      4.65 0.00000434
3 genderfemale   -0.198    0.0510     -3.88 0.000120  

The models we’re describing here contain only additions, not interactions. So, you may be wondering, “What if I want to model more complex relationships between predictors?” If that’s what you’re trying to do now, you may want to jump ahead to discovr_10 for moderation and mediation.

Otherwise, we will get there in this course eventually!

Comparing Models

Now we have two models, one simpler with only a single predictor, and the other with two predictors. We might next want to test whether the more complex, two-predictor model is a significant improvement over the simpler model, in order to decide which model to retain. We can do this with the anova() function1 to compare the two models.

Warning

The anova() function will only work for model comparison for particular models.

  1. The models must be hierarchical. That is, the more complex model(s) must contain all of the predictor(s) present in the less complex model(s).

  2. All models must be fit to the same dataset. If, for example, your first predictor has no missing values, but your second predictor had one, the model with using the second predictor would have been fit to a dataset of a different size than the model using only the first, and the anova() function will throw an error to this effect.

If you encounter this issue, you may need to inspect and clean your dataset before you proceed with analysis!

Exercise

Put both model objects into the anova() function to find out which model to retain.

anova(eval_lm, eval_full_lm)
Analysis of Variance Table

Model 1: eval ~ beauty
Model 2: eval ~ beauty + gender
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1    461 137.16                                  
2    460 132.81  1    4.3467 15.056 0.0001196 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The output from many base-R {stats} functions, like anova(), include a line labeled Signif. codes that provide a key for understanding the notation given for significance levels in p-values.

Reading the key from left to right, we can see that a result is given three asterisks (***) when the p-value is between 0 and .001; two asterisks between .001 and .01; and so on.

This can be a useful visual check, especially because p-values that are very, very small are frequently expressed in scientific notation, which can make them more difficult to spot.

The F-test is significant, indicating that the more complex two-predictor model is a significant improvement over the one-predictor model, so we will proceed with the two-predictor model.

Standardised Bs

You may have noticed that the output we’ve seen so for only contains unstandardised model parameter estimates. If we want standardised Bs expressing the relationship between predictor(s) and outcome in standard deviation units, we can make use of the model_parameters() function from the {parameters} package to standardise our bs.

Exercise

Use the standarize = "refit" argument in the parameters::model_parameters() function to obtain standardised Bs.

Warning

Note the spelling of standardize with a “z”! Spelling it with an “s” will not rescale the parameter estimates.

parameters::model_parameters(eval_full_lm, standardize = "refit")
Parameter       | Coefficient |   SE |         95% CI | t(460) |      p
-----------------------------------------------------------------------
(Intercept)     |        0.15 | 0.06 | [ 0.03,  0.27] |   2.53 | 0.012 
beauty          |        0.21 | 0.05 | [ 0.12,  0.30] |   4.65 | < .001
gender [female] |       -0.36 | 0.09 | [-0.54, -0.18] |  -3.88 | < .001

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation.

Assumptions Checks

At Sussex, we teach a range of model diagnostics and robust model sensitivity tests in order to test model assumptions. We will look briefly at each of these in turn.

Tip

Remember, there are more examples and longer explanations in the discovr_08 tutorial!

Residual Plots

To begin, we can generate nicely formatted, customisable residual plots using the function ggplot2::autoplot(). However, it is essential to load the {ggfortify} package in order for this to work!

Exercise

Load the {ggfortify} package and use the autoplot() function to generate residual plots for eval_full_lm. Set the which argument to c(1, 3, 2, 4).

library(ggfortify)
Warning: package 'ggfortify' was built under R version 4.3.1
Loading required package: ggplot2
Warning: package 'ggplot2' was built under R version 4.3.1
ggplot2::autoplot(eval_full_lm, which = c(1, 3, 2, 4))

If you’re wondering what’s up with which, the plot() function that autoplot() is based on has a total of six plots available. Here I’ve chosen the two residual plots, the normal Q-Q and the Cook’s distance plot. Plots 5 and 6 have to do with leverage, which we don’t teach at UG level. To see them, simply add them to the which = argument.

If you want to customise the theme or look of these plots further, they are built with {ggplot2} so you can add or change anything about them using that package. We will come round to a detailed exploration of data visualisations with {ggplot2} in Tutorial 06.

Distribution of Standardised Residuals

In the UG core statistics modules at Sussex, we teach that normality is the least important of the assumptions of the linear model - the most important being additivity and linearity - so we do not generally worry too much about normally distributed standardised residuals, especially in large sample sizes. However, we do teach them how to evaluate the proportion of standardised residuals with values above $$1.96 (approximately 5%), $$2.56 (approximately 1%), and above $$3 (likely an outlier). For details on how to use broom::augment() to obtain standardised residuals and other model diagnostic measures like Cook’s distance, see the discovr_08 tutorial.

Robust Models

At UG level, we teach robust models for two purposes:

  1. As sensitivity tests to check assumptions. If a robust technique that adjusts for a particular issue, such as heteroscedasticity, results in a model that is substantially different from the unadjusted model, we might conclude that the unadjusted model did in fact have that particular issue.
  2. As robust alternatives to the unadjusted model.
Robust Parameter Estimates

Our first robust model re-estimates the parameter estimates using robust techniques withe the robust::lmRob() function. We can then compare the robust parameter estimates to the unadjusted ones obtained with lm() to find out if our estimates were biased, and report the robust parameter estimates if they were.

Exercise

Fit the same two-predictor model again with robust::lmRob() and compare the results to the unadjusted two-predictor model.

eval_lm_rob <- robust::lmRob(eval ~ beauty + gender, data = teach_tib)

eval_lm_rob

Call:
robust::lmRob(formula = eval ~ beauty + gender, data = teach_tib)

Coefficients:
 (Intercept)        beauty  genderfemale  
      4.1044        0.1388       -0.2071  
eval_full_lm

Call:
lm(formula = eval ~ beauty + gender, data = teach_tib)

Coefficients:
 (Intercept)        beauty  genderfemale  
      4.0816        0.1486       -0.1978  

Comparing the values of the two versions of the model, we can see that the parameter estimates have changed very little, so we might conclude that the unadjusted model was fine.

Robust CIs and p-values

To test and adjust for heteroscedastic residuals, we can re-estimate the standard error using a robust method. To do this, we’ll use the parameters::model_parameters() function with the argument vcov = "HC4" as recommended in the discovr tutorial. To do this, use the unadjusted model object eval_full_lm as the first argument.

Exercise

Use parameters::model_parameters() to re-estimate the SEs, CIs, and p-values, and compare the results to the unadjusted model.

parameters::model_parameters(eval_full_lm, vcov = "HC4")
Parameter       | Coefficient |   SE |         95% CI | t(460) |      p
-----------------------------------------------------------------------
(Intercept)     |        4.08 | 0.03 | [ 4.02,  4.15] | 125.28 | < .001
beauty          |        0.15 | 0.03 | [ 0.08,  0.21] |   4.59 | < .001
gender [female] |       -0.20 | 0.05 | [-0.30, -0.10] |  -3.94 | < .001

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation.
broom::tidy(eval_full_lm, conf.int = TRUE)
# A tibble: 3 × 7
  term         estimate std.error statistic    p.value conf.low conf.high
  <chr>           <dbl>     <dbl>     <dbl>      <dbl>    <dbl>     <dbl>
1 (Intercept)     4.08     0.0329    124.   0            4.02      4.15  
2 beauty          0.149    0.0320      4.65 0.00000434   0.0858    0.211 
3 genderfemale   -0.198    0.0510     -3.88 0.000120    -0.298    -0.0976

The parameter estimates will not change, but the SEs, CIs, and p-values may. Here, the values are nearly identical and there are no major changes - that is, no predictors have become non-significant that were previously significant - so we might again conclude that the unadjusted model was not unduly biased.

Bootstrapping

If we had a small sample size, a final option would be to bootstrap the confidence intervals. To do this, we will again use parameters::model_parameters(), but this time with bootstrap = TRUE.

Note: Sample size is not an issue with this dataset (N = 463).

Exercise

Produce bootstrapped confidence intervals for the two-predictor model and compare to the unadjusted confidence intervals.

parameters::model_parameters(eval_full_lm, bootstrap = TRUE)
Parameter       | Coefficient |         95% CI |      p
-------------------------------------------------------
(Intercept)     |        4.08 | [ 4.02,  4.15] | < .001
beauty          |        0.15 | [ 0.08,  0.21] | < .001
gender [female] |       -0.20 | [-0.30, -0.10] | < .001

Uncertainty intervals (equal-tailed) are naıve bootstrap intervals.
broom::tidy(eval_full_lm, conf.int = TRUE)
# A tibble: 3 × 7
  term         estimate std.error statistic    p.value conf.low conf.high
  <chr>           <dbl>     <dbl>     <dbl>      <dbl>    <dbl>     <dbl>
1 (Intercept)     4.08     0.0329    124.   0            4.02      4.15  
2 beauty          0.149    0.0320      4.65 0.00000434   0.0858    0.211 
3 genderfemale   -0.198    0.0510     -3.88 0.000120    -0.298    -0.0976

We could once again note that there are no major changes, as previously, so the evidence of our checks suggests that the original, unadjusted model was not unduly biased.

Quarto

Quarto documents are a mix of text and code, the next generation of R Markdown documents. For our purposes, we will be using R within Quarto, but Quarto documents support the integration of many different coding languages, including Python, Julia, and Observable. If you’ve previously used Rmd (RMarkdown), Quarto is backwards-compatible and will able to render most documents with no issues.

Help with Quarto

The official Quarto Guide is extensive, detailed, and extremely helpful. It’s always the best first stop for any questions you have about using Quarto. Quarto also offers detailed tutorials.

This quick-reference to Markdown formatting is particularly helpful.

Prof Andy Field has also recorded a series of video guides to using Quarto that are used in UG teaching.

Getting Started

To get some hands-on practice working with Quarto, we will create a new Quarto document from scratch. We will use the analysis code we’ve already written to create a nicely formatted report.

Exercise

Create a new Quarto document via File > New File > Quarto Document.

If you like, you can give it a title; you will be able to change this later.

Then, click “Create”.

By default, your new Quarto document will already have some settings and content to demonstrate how it works. Most importantly, you can see that there are three main types of information in this document:

  1. The YAML header at the top, delineated by ---s, which contains information about how the document will be rendered
  2. The body text, which contains regular (i.e. non-code) text
  3. The code chunks, which contain code - in this case, specifically, R.

A screenshot of the default Quarto document, with the YAML, text, and code chunks each circled and labelled in a different colour.

Exercise

To complete our setup, take the following steps:

  • Delete everything in the new Quarto document except for the YAML header (i.e. all the text and code chunks).
  • Clear your Environment by clicking the broom icon in the Environment tab.
  • Restart your R session (via Session > Restart R).
Warning

Make sure you have completed and saved your worksheet with all of your code before you do this!

A new feature with Quarto is Visual mode, which is very like using Word or a similar word processing programme. All formatting can be applied and previewed in the document itself, using keyboard shortcuts or the familiar formatting buttons along the top toolbar.

Visual mode also has an “insert anything” shortcut, /, that allows you to quickly insert elements into your document. If you type / at the start of a new line (or Ctrl/Cmd + / otherwise), a drop-down list of possible elements will appear, which you can navigate by scrolling or search by typing.

Alternatively, you can toggle to Source mode by clicking the “Source” button at the very top left of the document pane. In Source mode you can edit the document using Markdown formatting. This allows more fine-tuned control of elements, but automatic formatting and the / shortcut don’t work in this mode.

Which you use is entirely personal preference - I toggle regularly between them depending on what I’m trying to do, so pick whichever works for you.

Creating a Code Chunk

To prepare for our analysis, we will need a place to load packages and read in the data. Any executable code - that is, code that you want to run and do something - must be written in a code chunk and NOT in the body text2.

Inserting Code Chunks

There are several ways to insert a new code chunk.

  • In Visual mode, by typing / and then Enter (since “R Code Chunk” is the first option)
  • In either mode, by:
    • Clicking the green “insert code chunk” button in the top right
    • Typing out the literal code fencing: ```{r} ```
    • Using a keyboard shortcut.
Exercise

Create a new code chunk and code that does the following:

  • Load all the packages we used in the previous Linear Model section
  • Read in the dataset located in the data folder of your project into an object called teach_tib
  • Create the two objects containing the linear models with one predictor and with two predictors.

Then, run the code chunk.

Load packages:

In the previous section, we used the following packages:

  • {tidyverse}, or specifically {broom}, for tidying up the linear model output
  • {ggfortify} for creating diagnostic plots
  • {parameters} for standardising bs, re-estimating CIs, and boostrapping
  • {robust} for robust parameter estimates

Read in the data:

In the Posit Cloud project workspace for this tutorial, there is a teaching.csv object. Look back on Tutorial 02 for more detail about how to read in datasets.

Create the models:

Copy and paste the code from your workbook document, or from the previous sections of this tutorial. You can also find all of the commands you have run your History tab (next to Environment).

Altogether, your new code chunk should look like this:

```{r}
library(tidyverse)
library(ggfortify)
library(parameters)
library(robust)

teach_tib <- readr::read_csv("data/teaching_ratings.csv")

eval_lm <- lm(eval ~ beauty, data = teach_tib)
eval_full_lm <- lm(eval ~ beauty + gender, data = teach_tib)
```

Run the code chunk by clicking the green “play” arrow, or by pressing Ctrl/Cmd + Shift + Enter while your cursor is inside the chunk.

Headings and Text

Headings

To begin, we’ll use headings to map out our document. Properly formatted headings are strongly recommended for any documents you write, for a variety of reasons:

  • They automatically create the outline (to the right) and navigation menu (to the bottom) of your document for easy navigation
  • They can be automatically converted into a table of contents
  • They are a crucial accessibility feature for navigating the document via keyboard/screenreader, as well as providing clear visual structure.
Inserting Headings

There are several ways to insert a new heading.

  • In Visual mode, by:
    • Using the text formatting dropdown. Click on “Normal” and select the heading level.
    • Using a keyboard shortcut, Ctrl + Alt + the number of the heading (e.g. 3 for a level 3 heading)
  • In either mode, by:
    • Typing hashes at the start of a new line, followed by a space (e.g. ### creates a level 3 heading)
Exercise

Create 2-3 level 2 headings in your document for the brief report of the linear model. You can choose anything you like, but the three I will refer to will be called “Model Comparison” (constructing and comparing the two models), “Assumptions Checks”, and “The Final Model” (reporting the final model in full).

Body Text

Any new-line text will automatically be plain text. This can be anything you like, although in our current document, we are writing a mock-formal results section.

Formatting Text

How to format body text depends on the mode you are in.

  • In Visual mode, by:
    • Using familiar keyboard shortcuts (e.g. Ctrl/Cmd + B for bold, Ctrl/Cmd + I for italics, etc.)
    • Using the formatting toolbar at the top of the document.
  • In Source mode, by using Markdown formatting.
Exercise

Under the first heading (which I have called “Model Comparison”), write a brief, journal-style description of the variables the two models contain.

If you aren’t inclined to write your own, here’s a brief sample text to use.

Two linear models were constructed to investigate the influences on teaching evaluation ratings. The first model contained only instructor beauty ratings as a predictor, and teaching evaluation ratings as the outcome (R2 = ???, F(???, ???) = ???, p = ???). The second model added instructor gender as a second predictor with no interaction, both again predicting teaching evaluation ratings (R2 = ???, F(???, ???) = ???, p = ???). An ANOVA comparing the models indicated a significant improvement in model fit for the second model compared to the first (F(???, ???) = ???, p = ???).

Dynamic Reporting

That last sentence should report which of the two models was better, based on the result of our F-test, rather than question marks. We could produce the output of this test, read it ourselves with our very own eyes/ears, and then type out the results by hand…but that is definitely not what we are going to do! Instead, we’ll look at a couple options for reporting the results dynamically. (We’ll come back to the model R2s in a moment.)

Inline Code

Our first option is to use inline code to report the numbers. In order to do this, we first need to know what information we have available in the test output, so we can make use of it.

Tip

This section will make extensive use of $ subsetting, which was covered in Tutorial 02, and [] subsetting, which was covered in Tutorial 01.

Exercise

In a new code chunk, use the broom::tidy() function to get a nice tibble of the anova() comparison between the two models, and save it in a new object called tidy_f. Then, print out this object in the Console.

We are creating the tidy_f object in a code chunk because we will want to use it again in our report.

We are calling this object in the Console because viewing its contents is only for our information/reference, and not something we want to appear or use directly in our finished document.

## In a code chunk
tidy_f <- broom::tidy(anova(eval_lm, eval_full_lm))

## In the Console
tidy_f
# A tibble: 2 × 7
  term                   df.residual   rss    df sumsq statistic   p.value
  <chr>                        <dbl> <dbl> <dbl> <dbl>     <dbl>     <dbl>
1 eval ~ beauty                  461  137.    NA NA         NA   NA       
2 eval ~ beauty + gender         460  133.     1  4.35      15.1  0.000120

The tibble we get here is quite different from the “Analysis of Variance Table” text output that we saw earlier when we printed out the results of the anova() function on its own. We now have all the values we need to report these results in a conveniently subsettable tibble with nice R-friendly names.

So, let’s have a go getting at some of those values.

Exercise

Get out the F-ratio of 15.0554899 from this object, and round it to two decimal places.

First, we need to get out the statistic variable, which we can do with $.

tidy_f$statistic
[1]       NA 15.05549

This returns a vector of two values, NA and 15.0554899. To index the second of these values, we need [] subsetting and the index number of the correct value.

tidy_f$statistic[2]
[1] 15.05549

Finally, we can use the round() function to round to two decimal places.

round(tidy_f$statistic[2], 2)
[1] 15.06

We now have a bit of R code that produces a single number that we want to report in the text3. The issue is that code chunks contain code, and body text contains text, but we would like the code that we’ve written to print out its value in the text!

The solution is inline code, a small bit of R code written in the body text (“inline”) that will, when the document is rendered, automatically insert the right number. Inline code is written as follows, in the text: `r some_r_code`

For our reporting, we will need to insert the inline code in exactly the spot where we would like the output of the code to appear. This particular bit of code produces the F-statistic, so we can replace the “???”s in our reporting where the value of the F-statistic should go:

An ANOVA comparing the models indicated a significant improvement in model fit for the second model compared to the first (F(???, ???) = `r round(tidy_f$statistic[2], 2)`, p = ???).

You can check whether this has worked in two ways:

  1. Place your cursor inside the backticks and press Ctrl/Cmd + Enter, as you would inside a code chunk. This will run the code and a small pop-up will show you what that code will produce when rendered.
  2. Render your document and see what appears!
Exercise

Use inline code to replace all of the ???s in the F-test reporting to produce a final report in APA style.

Hint: Rounding p-values is a bit tricky. Check out the {papaja} package to see if you can find a function besides round() that will produce the p-value with the correct formatting.

All of the other numbers should be straightforward, except for the p-value. Rounding to three decimal places with round() will result in a value of 0, which is not what we want. Instead, since the p-value is below .001, we want “< .001”.

{papaja} has the printp() function, which will do this exactly was we like (as well as containing a lot of other useful functions for rounding and printing in APA style!)

Your final text may look like this:

An ANOVA comparing the models indicated a significant improvement in model fit for the second model compared to the first (F(`r tidy_f$df[2]`, `r tidy_f$df.residual[2]`) = `r round(tidy_f$statistic[2], 2)`, p `r papaja::printp(tidy_f$p.value[2]`).

Note that there’s no need for a “<” symbol because papaja::printp() includes it automatically.

This will render as:

An ANOVA comparing the models indicated a significant improvement in model fit for the second model compared to the first (F(1, 460) = 15.06, p < .001).

CHALLENGE: Create a bit of inline code that will either report a significant or non-significant result depending on the value of p.

Hint: You may need to check out the ifelse() function.

If we really want our reporting to be resilient, we want to remove or replace all the places where we have to manually remember to update the code if our results change. In this case, our reporting reads:

An ANOVA comparing the models indicated a significant improvement in model fit…

But if we wanted to reuse this code for another report, we would have to remember to update this depending on the actual value of p. OR, we can have R do it for us.

First, we need to write a bit of R code that will evaluate whether the value of p is above or below a particular threshold, and then output the correct text. We could do this with ifelse(), a handy little base R function with three arguments. The first is a test returning a logical value (either TRUE or FALSE). If the test returns TRUE, the second argument is executed; if the test returns FALSE, the final argument is executed.

ifelse(
1  tidy_f$p.value[2] < .05,
2  "significant",
3  "non-significant"
)
1
The test: is the p-value for our F-test less than .05? (If you have a different \(\alpha\)-threshold, you could hard-code it here or use an object you’ve defined previously for this comparison.)
2
What to do if the test returns TRUE: print “significant”.
3
What to do if the test returns FALSE: print “non-significant”.
[1] "significant"

We can then take this command and replace the word “significant” in our report with this inline code:

An ANOVA comparing the models indicated a `r ifelse(tidy_f$p.value[2] < .05, "significant", "non-significant")` improvement in model fit…

Now we don’t have to worry about getting this right: our document, when rendered, will automatically insert the right word depending on the data.

Automatic Reporting

As helpful as inline code is (and I would recommend reporting all values dynamically/automatically wherever possible, so it is very useful!), you may have noticed that there was a lot of repetitive typing that also made the text itself quite difficult to read, as well lots of opportunities make typos or mistakes. Surely there’s a simpler way to do this sort of thing!

There are, in fact, many simpler ways to do common tasks like this, which take advantage of the fact that an object created by a particular function will always have the same structure. One option is to make further use of the {papaja} package, which is designed for just this purpose.

Exercise

Use the {papaja} documentation to fill in the statistical reporting for each of the linear models (i.e., R2, F, and p) using only one piece of inline code for each.

When you’re done, render your document to see the results!

The {papaja} documentation illustrates the process using t.test(), which works the same way as lm(). The key here is to use the original objects containing the models you want to report.

Let’s have a look at the first of the two models and see what the papaja::apa_print() function gives us.

papaja::apa_print(eval_lm)
$estimate
$estimate$Intercept
[1] "$b = 4.00$, 95\\% CI $[3.95, 4.05]$"

$estimate$beauty
[1] "$b = 0.13$, 95\\% CI $[0.07, 0.20]$"

$estimate$modelfit
$estimate$modelfit$r2
[1] "$R^2 = .04$"

$estimate$modelfit$r2_adj
[1] "$R^2_{adj} = .03$"

$estimate$modelfit$aic
[1] "$\\mathrm{AIC} = 756.65$"

$estimate$modelfit$bic
[1] "$\\mathrm{BIC} = 769.06$"



$statistic
$statistic$Intercept
[1] "$t(461) = 157.73$, $p < .001$"

$statistic$beauty
[1] "$t(461) = 4.13$, $p < .001$"

$statistic$modelfit
$statistic$modelfit$r2
[1] "$F(1, 461) = 17.08$, $p < .001$"



$full_result
$full_result$Intercept
[1] "$b = 4.00$, 95\\% CI $[3.95, 4.05]$, $t(461) = 157.73$, $p < .001$"

$full_result$beauty
[1] "$b = 0.13$, 95\\% CI $[0.07, 0.20]$, $t(461) = 4.13$, $p < .001$"

$full_result$modelfit
$full_result$modelfit$r2
[1] "$R^2 = .04$, $F(1, 461) = 17.08$, $p < .001$"



$table
A data.frame with 6 labelled columns:

       term estimate     conf.int statistic  df p.value
1 Intercept     4.00 [3.95, 4.05]    157.73 461  < .001
2    Beauty     0.13 [0.07, 0.20]      4.13 461  < .001

term     : Predictor 
estimate : $b$ 
conf.int : 95\\% CI 
statistic: $t$ 
df       : $\\mathit{df}$ 
p.value  : $p$ 
attr(,"class")
[1] "apa_results" "list"       

We’ve got a huge number of options here, but for this exercise we wanted R2, F, and p. All three are given under $full_result$modelfit$r2. We will need to save the output from apa_print() into an object, then we can subset it using inline code:

eval_lm_out <- papaja::apa_print(eval_lm)
eval_full_lm_out <- papaja::apa_print(eval_full_lm)

Two linear models were constructed to investigate the influences on teaching evaluation ratings. The first model contained only instructor beauty ratings as a predictor, and teaching evaluation ratings as the outcome (`r eval_lm_out$full_result$modelfit$r2`). The second model added instructor gender as a second predictor with no interaction, both again predicting teaching evaluation ratings (`r eval_full_lm_out$full_result$modelfit$r2`).

Which will render as:

Two linear models were constructed to investigate the influences on teaching evaluation ratings. The first model contained only instructor beauty ratings as a predictor, and teaching evaluation ratings as the outcome (\(R^2 = .04\), \(F(1, 461) = 17.08\), \(p < .001\)). The second model added instructor gender as a second predictor with no interaction, both again predicting teaching evaluation ratings (\(R^2 = .07\), \(F(2, 460) = 16.33\), \(p < .001\)).

Table Formatting

Next, we will jump to the “Final Model” heading and have a look at how to turn our final model output into a nicely formatted table. Once again, {papaja} provides a quick and beautiful solution for reporting, so let’s use it again.

Exercise

Using the {papaja} help documentation, produce a nicely formatted table of the final model, presenting the parameter estimates, p-values etc. for each predictor under the third (“Final Model”) heading.

We already have the necessary object, eval_full_lm_out, from the previous task. We just need to subset it as described in the help documentation.

This command should go in a new code chunk, wherever you want the table to appear in your document.

papaja::apa_table(eval_full_lm_out$table)
(#tab:unnamed-chunk-19)
Predictor \(b\) 95% CI \(t\) \(\mathit{df}\) \(p\)
Intercept 4.08 [4.02, 4.15] 123.94 460 < .001
Beauty 0.15 [0.09, 0.21] 4.65 460 < .001
Genderfemale -0.20 [-0.30, -0.10] -3.88 460 < .001

CHALLENGE: {papaja} isn’t the only package to provide easy formatting for commonly reported tests. Have a go creating this table again using the nice_table() function from the {rempsyc} package, which allows a bit more flexibility in

The nice_table() function can be for tables generally, but it can apply specialised formatting for model tables created with broom, if we use the broom = argument to tell the function what formatting template to apply.

rempsyc::nice_table(broom::tidy(eval_full_lm, conf.int = TRUE), broom = "lm")

Term

b

SE

t

p

95% CI

(Intercept)

4.08

0.03

123.94

< .001***

[4.02, 4.15]

beauty

0.15

0.03

4.65

< .001***

[0.09, 0.21]

genderfemale

-0.20

0.05

-3.88

< .001***

[-0.30, -0.10]

Cross-Referencing

As anyone who has had to create a long document with lots of tables and figures knows, keeping track of the numbering is a huge pain, especially when, for instance, a reviewer asks you to add something partway through and then everything has to be renumbered4.

The good news is that Quarto can take care of figure and table numbering automatically. There are two steps to this:

  • Include a label in the relevant code chunk, using the prefix fig- for figures and tbl- for tables.
  • Refer to the figure or table in the text using @.
Exercises

Using the Quarto help documentation, write a short introductory sentence under the “Final Model” heading and refer to the final model table with a cross-reference.

First, add a label and caption to the code chunk from the previous task that produces the model table.

```{r}
#| label: tbl-final-model
#| tbl-cap: "The final model predicting teaching evaluation ratings from instructor beauty and gender."

papaja::apa_table(eval_full_lm_out$table,
                  caption = "The final model predicting teaching evaluation ratings from instructor beauty and gender.")
```

You can of course write whatever you like, or borrow the text below, but use “Table 1” (or whatever label your gave the table code chunk) to refer to the table.

The final model with two predictors is presented in full in @tbl-final-model.

Altogether, it should render as follows:

The Final Model

The final model with two predictors is presented in full in Table 1.

papaja::apa_table(eval_full_lm_out$table)
Table 1: The final model predicting teaching evaluation ratings from instructor beauty and gender.
Predictor \(b\) 95% CI \(t\) \(\mathit{df}\) \(p\)
Intercept 4.08 [4.02, 4.15] 123.94 460 < .001
Beauty 0.15 [0.09, 0.21] 4.65 460 < .001
Genderfemale -0.20 [-0.30, -0.10] -3.88 460 < .001

CHALLENGE: Complete the final “Assumptions Checks” section summarising the checks and using figure cross-referencing to insert and refer to the diagnostic plots.

Hint: To report the exact maximum value of Cook’s distance, you will also need to refer to discovr_08 for how to use broom::augment().

You can write whatever you like, but here’s a suggestion with the plot included.

Assumptions Checks

We next assessed the model with two predictors for any evidence of bias. Residual plots (@fig-diag-plots) did not indicate any outstanding issues with normality, linearity, or heteroscedasticity. There was also no evidence of influential cases, as the max value of Cook’s distance was `r round(max(broom::augment(eval_full_lm)$.cooksd), 2)`. Robust models were also fitted as sensitivity checks. Robust parameter estimates estimated using the {robust} package were minimally different from the unadjusted parameter estimates. Similarly, robust HC4 standard errors estimated using the {parameters} package yielded confidence intervals and p-values very similar to the unadjusted values. Therefore, we will proceed with the unadjusted two-predictor model as our final model.

```{r}
#| label: fig-diag-plots
#| fig-cap: "Diagnostic plots for the two-predictor model."

ggplot2::autoplot(eval_full_lm, which = c(1, 3, 2, 4))
```

Rendering

After all this work to analyse, interpret, and report, it’s finally time to produce the final document. The process of turning a Quarto document into some output format - including running all the code and applying all the formatting - is called rendering (previously knitting with RMarkdown).

Exercise

Render your report document using the “Render” button at the top of the document, or by using the keyboard shortcut Ctrl/Cmd + Shift + K5.

Global Options

At the moment, your report may not be as clean as we’d like it to be: there are likely messages from R floating around and code all over the place, whereas for a formal report, we of course only want to show the final output. This is where the YAML header comes in.

The default YAML header contains only a few explicit settings, and if you haven’t changed anything, probably looks like this:

---
title: "Untitled"
format: html
editor: visual
---

By adding options to the YAML header, we can determine how the document as a whole is rendered. Here are the common ones that I use on the regular:

---
title: "Linear Model Report"
format:
  html:
1    toc: true
editor: visual
2self-contained: true
3execute:
4  echo: false
5  warning: false
6  message: false
---
1
Automatically produce a table of contents (ToC) from the document headings.
2
Combine all the files, images, stylesheets etc. into a single output document. Necessary if you want to send an HTML file to someone else and have it look as it should!
3
Set default behaviour for all code chunks as follows
4
Run code and show output, but do not show the code itself.
5
Do not show any warnings produced by code.
6
Do not show any messages produced by code.

The requirements for each document will change depending on its purpose, so you will likely find that

Tip

The Quarto help documentation, as usual, has a complete list of YAML options, including how to set default behaviour for figures and other settings.

Code Chunk Options

Global options apply to the entire document, but you may want to change these settings for individual code chunks to override the default settings in the YAML.

For example, you may have a code chunk containing some processing code that you used to view and clean your data. If your global execution option is set to echo: false, the code output would still appear, although the code itself would be hidden. If, for this particular code chunk, you don’t want the output to appear either, you can override the global option with a local option for that chunk only.

Local code chunk options appear as “hashpipe” (#|) comments, as we have seen earlier with labels. They use the same syntax as YAML, but the settings only apply to individual code chunks. Any settings that aren’t explicitly changed within the chunk are inherited from the YAML settings.

For this example, we could set a local code chunk option include: false which will prevent the output from appearing in the document.6

```{r}
#| include: false

some_code_doing_cleaning_and_processing
```

Output Formats

For these tutorials, we will generally stick to HTML, as it’s the most painless of the rendering options. However, you will likely find yourself wanting to produce some other type of document, which you can easily7 do from the same Quarto document.

To render to a different format, change the YAML format: setting to a different output.

Tip

As per, the Quarto guide on output formats has all the information you need!

Exercise

Render your linear model report to a Word document.

Simply update the format: html YAML option to format: docx and render your document. Note that the format options are usually named after the file extension rather than the name of the programme necessarily.

Well done!

And there we have it! You now have a complete example of a linear model report, rendered into both HTML and Word, to refer to. It’s amazing how much you’re able to do after just a few weeks!

This is the end of the FundRmentals section of the course. If you’re so inclined, we’ll see you in the Essentials section, which will cover data wrangling and cleaning, and running and reporting many more statistical analyses.

Footnotes

  1. Somewhat confusingly, this is not the function we teach UGs to perform ANOVAs! See discovr tutorials 11, 12, 15, and 16 for a detailed guide through the afex package for running linear models with categorical predictors.↩︎

  2. Mostly true, except for inline code!↩︎

  3. Generally, inline code should only ever produce a single value, otherwise the formatting can get interesting. This value, however, could longer than a single number, if you want to get creative with your dynamic reporting!↩︎

  4. ↩︎

  5. It seems obvious that this should be R instead of K, but remember this made perfect sense when it was called “knitting”!↩︎

  6. For writing tutorials, I make extensive use of eval: false, which includes the code in the output but does not attempt to run the code. This allows me to write all kinds of nonsense code without R getting stroppy and throwing errors all of the place!↩︎

  7. Depends crucially on what you consider to be “easy”, especially when dealing with PDFs!↩︎